library(plyr)
library(tidyverse)
library(magrittr)
library(emmeans)
library(stringr)
library(broom)
library(stargazer)
library(gtools)
library(purrrlyr)
library(ggstance)

# first set wd variable to location of d1, d2, & d3
# wd <- "C:/Dropbox/Babak Ethan Tom/data for r&p/"

d1 <- str_c(wd, "d1.RDS") %>% 
  readRDS


d2 <- str_c(wd, "d2.RDS") %>% 
  readRDS

d3 <- str_c(wd, "d3.RDS") %>% 
  readRDS

d4 <- str_c(wd, "d4.RDS") %>% 
  readRDS

lm_3 <- map2(
  str_c(
    c("polar_num",
      "reg_num",
      "paris_num",
      "reg_num"),
    " ~ cond"
  ),
  list(
    d1, d1, d2, d2
  ),
  function(i, j)
    lm(i, j)
)

f1 <- pmap_dfr(
  list(
    lm_3,
    list(
      d1, d1, d2, d2
    ),
    c("polar_num",
      "reg_num",
      "paris_num",
      "reg_num")
  ),
  function(i, j, k)
    # i <- lm_3[[2]]
    # j <- d1
    # k <- "reg_num"
    
    i %>%
    emmeans(
      revpairwise ~ cond, 
      data = j
    ) %>% 
    extract2(1) %>% 
    tidy %>%
    mutate(
      dv = k
    )
) %>% 
  mutate(
    misstatement = c("Trump claims polar ice-caps at all time high",
                     "Trump claims Paris accords impose requirements only on US") %>% 
      rep(, each = 6),
    dv = dv %>% 
      mapvalues(
        c("paris_num",
          "polar_num",
          "reg_num"),
        c("Paris Accords allow all states flexibility in meeting targets (5pt-scale--larger values indicate accuracy)",
          "Polar ice caps at all time low  (5pt-scale--larger values indicate accuracy)",
          "Favor regulations to protect environment (3pt-scale--larger values indicate support for regulations)")
      ) %>% 
      factor(
        c("Paris Accords allow all states flexibility in meeting targets (5pt-scale--larger values indicate accuracy)",
          "Polar ice caps at all time low  (5pt-scale--larger values indicate accuracy)",
          "Favor regulations to protect environment (3pt-scale--larger values indicate support for regulations)")
      ),
    lab = estimate %>% 
      round(2)
  )

f1_l <- pmap_dfr(
  list(
    lm_3,
    list(
      d1, d1, d2, d2
    ),
    c("polar_num",
      "reg_num",
      "paris_num",
      "reg_num")
  ),
  function(i, j, k)
    # i <- lm_3[[2]]
    # j <- d1
    # k <- "reg_num"
    
    i %>%
    emmeans(
      revpairwise ~ cond, 
      data = j
    ) %>% 
    extract2(2) %>% 
    tidy %>%
    mutate(
      dv = k
    )
) %>% 
  filter(
    p.value < .05
  ) %>% 
  mutate(
    lab = str_c(
      "p=",
      p.value %>% 
        signif(digits = 1) %>% 
        format(scientific = F) %>% 
        as.character %>% 
        str_sub(2)
    )
  )

f1_l$lab %<>% 
  str_sub(
    , 
    f1_l$lab %>% 
      map_int(
        function(i)
          i %>% 
          str_locate(
            1:9 %>% 
              str_c(collapse = "|")
          ) %>% 
          extract(1)
      )
  )


# figure 1
p1 <- f1 %>% 
  ggplot() +
  geom_pointrange(
    aes(misstatement, estimate, color = cond,
        group = cond,
        ymin = conf.low, 
        ymax = conf.high),
    position = position_dodge(width = .6),
    shape = 21,
    size = 1.1,
    fatten = 6.5,
    fill = "white"
  ) +
  geom_text(
    aes(misstatement, estimate, label = lab,
        group = cond),
    position = position_dodge(width = .6),
    size = 2.5
  ) +
  geom_segment(
    aes(x = 1, xend = 1.2,
        y = 3.3, yend = 3.3),
    linetype = "solid",
    size = .5,
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(2)
    )
  ) +
  geom_segment(
    aes(x = 1.2, xend = 1.2,
        y = 3.3, yend = 3.71),
    linetype = "solid",
    size = .5,
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(2)
    )
  ) +
  geom_segment(
    aes(x = 1, xend = 1,
        y = 3.3, yend = 3.36),
    linetype = "solid",
    size = .5,
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(2)
    )
  ) +
  geom_text(
    aes(x = 1.1, 
        y = 3.25),
    label = f1_l$lab[[2]],
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(2)
    ),
    size = 3
    ) +
  geom_segment(
    aes(x = 1, xend = .8,
        y = 3.1, yend = 3.1),
    linetype = "solid",
    size = .5,
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(2)
    )
  ) +
  geom_segment(
    aes(x = 1, xend = 1,
        y = 3.1, yend = 3.2),
    linetype = "solid",
    size = .5,
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(2)
    )
  ) +
  geom_segment(
    aes(x = .8, xend = .8,
        y = 3.1, yend = 3.59),
    linetype = "solid",
    size = .5,
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(2)
    )
  ) +
  geom_text(
    aes(x = .9, 
        y = 3.05),
    label = f1_l$lab[[1]],
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(2)
    ),
    size = 3
  ) +
  geom_segment(
    aes(x = 1, xend = 1.2,
        y = 2.825, yend = 2.825),
    linetype = "solid",
    size = .5,
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(1)
    )
  ) +
  geom_segment(
    aes(x = 1.2, xend = 1.2,
        y = 2.825, yend = 3.21),
    linetype = "solid",
    size = .5,
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(1)
    )
  ) +
  geom_segment(
    aes(x = 1, xend = 1,
        y = 2.825, yend = 2.92),
    linetype = "solid",
    size = .5,
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(1)
    )
  ) +
  geom_text(
    aes(x = 1.1, 
        y = 2.825 -.05),
    label = f1_l$lab[[3]],
    data = data.frame(
      dv = f1$dv %>% 
        levels %>% 
        extract(1)
    ),
    size = 3
  ) +
  scale_color_grey(
    start = .01,
    end = .8,
    guide = guide_legend(
      title = "",
      override.aes = list(
        shape = 16
      )
    ),
    breaks = f1$cond %>% 
      levels,
    labels= c(
      "Survey items only",
      "Items &\nTrump Misstatement",
      "Items, Misstatement,\n& Correction"
    ) 
  ) +
  scale_x_discrete(
    breaks =  f1$misstatement %>%
      factor %>% 
      levels,
    labels= f1$misstatement %>%
      factor %>% 
      levels %>% 
      str_wrap(width = 25)
  ) +
  facet_grid(
    .~dv, 
    scales = "free_x",
    space = "free", 
    labeller = label_wrap_gen(30)
  ) +
  labs(
    x = "",
    y = "",
    colors = "",
    fill = "",
    group = ""
  ) +
  theme_minimal() +
  theme(
    panel.background  =  element_rect(
      color =  "grey95",
      fill = "grey95"
    ),
    panel.grid = element_blank(),
    strip.background = element_rect(
      color =  "grey95",
      fill = "grey95"
    ),
    axis.text.y = element_text(hjust = .5),
    strip.text.y = element_text(angle = 0),
    legend.background = element_rect(color = "white",
                                     fill = "white"),
    legend.margin = margin(-.4, 0, 0, 0, "cm"),
    panel.grid.minor = element_blank(),
    plot.title = element_text(face = "bold", hjust = 0),
    plot.caption = element_text(face = "italic"),
    legend.position = "bottom"
    ) 


# table 2

lm_3 %>% 
  stargazer(
    add.lines = list(
      "Auxiliary Quantitities",
      c("Misstatement - Correction",
        map2_chr(
          lm_3,
          list(
            d1, d1, d2, d2
          ),
          function(i, j){
            # i <- lm_3[[1]]
            # j <- d1
            
            k <- i %>%
              emmeans(pairwise ~ cond, data = j) %>% 
              use_series(contrasts) %>%
              as.data.frame %>% 
              filter(
                contrast == "misstatement - correction"
              )
            
            str_c(
              k$estimate %>% 
                round(2) %>% 
                as.character %>% 
                str_replace("0.", "."),
              k$p.value %>% 
                stars.pval(),
              " (",
              k$SE %>% 
                round(2) %>% 
                as.character %>% 
                str_replace("0.", "."),
              ")"
            )
          }
        )
      ),
      c("Items - Correction",
        map2_chr(
          lm_3,
          list(
            d1, d1, d2, d2
          ),
          function(i, j){
            # i <- lm_3[[1]]
            # j <- d1
            
            k <- i %>%
              emmeans(pairwise ~ cond, data = j) %>% 
              use_series(contrasts) %>%
              as.data.frame %>% 
              filter(
                contrast == "items - correction"
              )
            
            str_c(
              k$estimate %>% 
                round(2) %>% 
                as.character %>% 
                str_replace("0.", "."),
              k$p.value %>% 
                stars.pval(),
              " (",
              k$SE %>% 
                round(2) %>% 
                as.character %>% 
                str_replace("0.", "."),
              ")"
            )
          }
        ) 
      ),
      c("Items - Misstatement",
        map2_chr(
          lm_3,
          list(
            d1, d1, d2, d2
          ),
          function(i, j){
            # i <- lm_3[[1]]
            # j <- d1
            
            k <- i %>%
              emmeans(pairwise ~ cond, data = j) %>% 
              use_series(contrasts) %>%
              as.data.frame %>% 
              filter(
                contrast == "items - misstatement"
              )
            
            str_c(
              k$estimate %>% 
                round(2) %>% 
                as.character %>% 
                str_replace("0.", "."),
              k$p.value %>% 
                stars.pval(),
              " (",
              k$SE %>% 
                round(2) %>% 
                as.character %>% 
                str_replace("0.", "."),
              ")"
            )
          }
        )
      )
    ),
    covariate.labels = c(
      "Condition: Misstatment",
      "Condition: Correction",
      "Constant"
    ),
    type = "text",
    # type = "latex",
    no.space = T,
    model.numbers = F,
    initial.zero = F,
    keep.stat = c("n", "rsq"),
    digits = 2,
    column.labels = c(
      "May Experiment",
      "June Experiment"
    ),
    column.separate = c(2, 2), 
    dep.var.labels.include = T,
    dep.var.labels   = 
      c("Polar caps size",
        "Enviro regs.",
        "Paris accords flex",
        "Enviro regs")
  )

# regression models for partisanship effects 

lm_1 <- str_c(
  c("polar_num",
    "reg_num"),
  " ~ cond * party_num"
  ) %>% 
  map(
    function(i)
      lm(
        i, data = d1
      )
  )

cont_1 <- map2_dfr(
  lm_1,
  c("polar caps",
    "regulations"),
  function(i, j)
    i %>% 
    emmeans(
      revpairwise ~ cond | party_num,
      at = list(
        party_num = seq(1, 5, length.out = 100)
      )
    ) %>% 
    extract2(2) %>% 
    tidy %>% 
    mutate(
      dv = j
    )
  ) %>% 
  mutate(
    cont2 = str_c(
      level1, "_", level2
    ) %>% 
      mapvalues(
        str_c(
          c("misstatement",
            "correction",
            "correction"),
          "_",
          c("items",
            "items",
            "misstatement")
        ),
        c("Misinformation effect",
          "Correction effect",
          "Difference in effects")
      ) %>% 
      factor(
        c("Misinformation effect",
          "Correction effect",
          "Difference in effects")
      ),
    lo = estimate %>% 
      subtract(
        std.error %>% 
          multiply_by(1.96)
      ),
    hi = estimate %>% 
      add(
        std.error %>%
          multiply_by(1.96)
      ),
    dv = dv %>% 
      mapvalues(
        c("polar caps",
          "regulations"),
        c("Polar ice caps' size at all time low (5pt scale)",
          "Favor regulation to protect environment (3pt scale)") %>% 
          str_wrap(width = 15)
      ) %>%
      factor(
        c("Polar ice caps' size at all time low (5pt scale)",
          "Favor regulation to protect environment (3pt scale)") %>% 
          str_wrap(width = 15)
      )
  )

c1b <- cont_1 %>% 
  group_by(dv, cont2) %>% 
  slice(c(1, 50, 100)) %>% 
  mutate(
    lab = estimate %>% 
      round(2) %>% 
      str_sub(2) %>% 
      as.character %>% 
      str_replace(fixed("0."), ".") %>% 
      str_pad(width = 3, side = "right", pad = "0"),
    lab = estimate %>% 
      sign %>% 
      equals(-1) %>% 
      ifelse(
        str_c(
          "-", lab
        ),
        lab
      ) %>% 
      str_c(
        c("***", "**", "*", "") %>% 
          extract(
            p.value %>% 
              findInterval(
                c(-Inf, .001, .01, .05, Inf)
              )
          )
      ),
    party_num = c(1, 3, 5)
  ) 

# figure 2
cont_1 %>% 
  ggplot() +
  geom_hline(
    yintercept = 0,
    linetype = "dashed"
  )  +
  geom_ribbon(
    aes(party_num,
        ymin = lo,
        ymax = hi),
    alpha = .4,
    size = .2,
    color = "grey5",
    fill = "grey90"
  ) +
  geom_segment(
    aes(x =  .6, xend = .6,
        y = .02 , yend = .75),
    arrow = arrow(
      type = "closed",
      length = unit(0.25, "cm")
    ),
    data = data.frame(
      dv = cont_1$dv %>% 
        levels %>% 
        extract2(1),
      cont2 = cont_1$cont2 %>% 
        levels %>% 
        extract2(1)
    ),
    color = "grey25"
  ) +
  geom_segment(
    aes(x =  .6, xend = .6,
        y = -.02 , yend = -.75),
    arrow = arrow(
      type = "closed",
      length = unit(0.25, "cm")
    ),
    data = data.frame(
      dv = cont_1$dv %>% 
        levels %>% 
        extract2(1),
      cont2 = cont_1$cont2 %>% 
        levels %>% 
        extract2(1)
    ),
    color = "grey25"
  ) +
  geom_segment(
    aes(x =  .6, xend = .6,
        y = -.02 , yend = -.25),
    arrow = arrow(
      type = "closed",
      length = unit(0.25, "cm")
    ),
    data = data.frame(
      dv = cont_1$dv %>% 
        levels %>% 
        extract2(2),
      cont2 = cont_1$cont2 %>% 
        levels %>% 
        extract2(1)
    ),
    color = "grey25"
  ) +
  geom_segment(
    aes(x =  .6, xend = .6,
        y = .02 , yend = .3),
    arrow = arrow(
      type = "closed",
      length = unit(0.25, "cm")
    ),
    data = data.frame(
      dv = cont_1$dv %>% 
        levels %>% 
        extract2(2),
      cont2 = cont_1$cont2 %>% 
        levels %>% 
        extract2(1)
      ),
    color = "grey25"
  ) +
  geom_text(
    aes(x, y, label = lab),
    size = 2.5,
    hjust = 0,
    angle = 90, 
    lineheight = .75,
    data = data.frame(
      dv = cont_1$dv %>% 
        levels %>% 
        extract2(1),
      cont2 = cont_1$cont2 %>% 
        levels %>% 
        extract2(1),
      x = .35,
      y = .05,
      lab = c("Respondents\nmore accurate")
      )
    ) +
  geom_text(
    aes(x, y, label = lab),
    size = 2.5,
    hjust = 0,
    angle = 90,
    lineheight = .75,
    data = data.frame(
      dv = cont_1$dv %>% 
        levels %>% 
        extract2(1),
      cont2 = cont_1$cont2 %>% 
        levels %>% 
        extract2(1),
      x = .35,
      y = -.6,
      lab = c("Respondents\nless accurate")
      )
    ) +
  geom_point(
    aes(party_num, estimate),
    data = c1b,
    shape = 21,
    size = 8,
    fill = "grey99",
    color = "grey5"
  ) +
  geom_text(
    aes(party_num, estimate, label = lab),
    data = c1b,
    size = 2.5
    ) +
  geom_text(
    aes(x, y, label = lab),
    size = 2.5,
    hjust = 0,
    angle = 90, 
    lineheight = .75,
    data = data.frame(
      dv = cont_1$dv %>% 
        levels %>% 
        extract2(2),
      cont2 = cont_1$cont2 %>% 
        levels %>% 
        extract2(1),
      x = .35,
      y = .025,
      lab = c("Respondents\nfavor regulation")
      )
    ) +
  geom_text(
    aes(x, y, label = lab),
    size = 2.5,
    hjust = 0,
    angle = 90,
    lineheight = .75,
    data = data.frame(
      dv = cont_1$dv %>% 
        levels %>% 
        extract2(2),
      cont2 = cont_1$cont2 %>% 
        levels %>% 
        extract2(1),
      x = .35,
      y = -.2,
      lab = c("Respondents\noppose regs.")
      )
    ) +
  facet_grid(
    dv ~ cont2,
    scales = "free",
    space = "free_x"
  ) +
  scale_x_continuous(
    expand = c(.08, .08),
    breaks = c(1, 3, 5),
    labels = c("Strong\nDem",
               "Independent",
               "Strong\nRep")
  ) +
  labs(
    x = "Partisanship",
    y = ""
  )

# figure 3

lm_2 <- str_c(
  c("paris_num",
    "reg_num"),
  " ~ cond * party_num"
  ) %>% 
  map(
    function(i)
      lm(
        i, data = d2
      )
  )


cont_2 <- map2_dfr(
  lm_2,
  c("paris coal",
    "regulations"),
  function(i, j)
    i %>% 
    emmeans(
      revpairwise ~ cond | party_num,
      at = list(
        party_num = seq(1, 5, length.out = 100)
      )
    ) %>%
    extract2(2) %>% 
    tidy %>% 
    mutate(
      dv = j
    )
) %>% 
  mutate(
    cont2 = str_c(
      level1, "_", level2
    ) %>% 
      mapvalues(
        str_c(
          c("misstatement",
            "correction",
            "correction"),
          "_",
          c("items",
            "items",
            "misstatement")
        ),
        c("Misinformation effect",
          "Correction effect",
          "Difference in effects")
      ) %>% 
      factor(
        c("Misinformation effect",
          "Correction effect",
          "Difference in effects")
      ),
    lo = estimate %>% 
      subtract(
        std.error %>% 
          multiply_by(1.96)
      ),
    hi = estimate %>% 
      add(
        std.error %>%
          multiply_by(1.96)
      ),
    dv = dv %>% 
      mapvalues(
        c("paris coal",
          "regulations"),
        c("Paris accords gave US flexibility on meeting targets (5pt scale)",
          "Favor regulation to protect environment (3pt scale)") %>% 
          str_wrap(width = 15)
      ) %>%
      factor(
        c("Paris accords gave US flexibility on meeting targets (5pt scale)",
          "Favor regulation to protect environment (3pt scale)") %>% 
          str_wrap(width = 15)
      )
  )

c2b <- cont_2 %>% 
  group_by(dv, cont2) %>% 
  slice(c(1, 50, 100)) %>% 
  mutate(
    lab = estimate %>% 
      round(2) %>% 
      str_sub(2) %>% 
      as.character %>% 
      str_replace(fixed("0."), ".") %>% 
      str_pad(width = 3, side = "right", pad = "0"),
    lab = estimate %>% 
      sign %>% 
      equals(-1) %>% 
      ifelse(
        str_c(
          "-", lab
        ),
        lab
      ) %>% 
      str_c(
        c("***", "**", "*", "") %>% 
          extract(
            p.value %>% 
              findInterval(
                c(-Inf, .001, .01, .05, Inf)
              )
          )
      ),
    party_num = c(1, 3, 5)
  ) 


cont_2 %>% 
  ggplot() +
  geom_hline(
    yintercept = 0,
    linetype = "dashed"
  )  +
  geom_ribbon(
    aes(party_num,
        ymin = lo,
        ymax = hi),
    alpha = .4,
    size = .2,
    color = "grey5",
    fill = "grey90"
  ) +
  geom_segment(
    aes(x =  .6, xend = .6,
        y = .02 , yend = .75),
    arrow = arrow(
      type = "closed",
      length = unit(0.25, "cm")
    ),
    data = data.frame(
      dv = cont_2$dv %>% 
        levels %>% 
        extract2(1),
      cont2 = cont_2$cont2 %>% 
        levels %>% 
        extract2(1)
    ),
    color = "grey25"
  ) +
  geom_segment(
    aes(x =  .6, xend = .6,
        y = -.02 , yend = -.75),
    arrow = arrow(
      type = "closed",
      length = unit(0.25, "cm")
    ),
    data = data.frame(
      dv = cont_2$dv %>% 
        levels %>% 
        extract2(1),
      cont2 = cont_2$cont2 %>% 
        levels %>% 
        extract2(1)
    ),
    color = "grey25"
  ) +
  geom_segment(
    aes(x =  .6, xend = .6,
        y = -.02 , yend = -.25),
    arrow = arrow(
      type = "closed",
      length = unit(0.25, "cm")
    ),
    data = data.frame(
      dv = cont_2$dv %>% 
        levels %>% 
        extract2(2),
      cont2 = cont_2$cont2 %>% 
        levels %>% 
        extract2(1)
    ),
    color = "grey25"
  ) +
  geom_segment(
    aes(x =  .6, xend = .6,
        y = .02 , yend = .3),
    arrow = arrow(
      type = "closed",
      length = unit(0.25, "cm")
    ),
    data = data.frame(
      dv = cont_2$dv %>% 
        levels %>% 
        extract2(2),
      cont2 = cont_2$cont2 %>% 
        levels %>% 
        extract2(1)
    ),
    color = "grey25"
  ) +
  geom_text(
    aes(x, y, label = lab),
    size = 2.5,
    hjust = 0,
    angle = 90, 
    lineheight = .75,
    data = data.frame(
      dv = cont_2$dv %>% 
        levels %>% 
        extract2(1),
      cont2 = cont_2$cont2 %>% 
        levels %>% 
        extract2(1),
      x = .35,
      y = .05,
      lab = c("Respondents\nmore accurate")
    )
    ) +
  geom_text(
    aes(x, y, label = lab),
    size = 2.5,
    hjust = 0,
    angle = 90,
    lineheight = .75,
    data = data.frame(
      dv = cont_2$dv %>% 
        levels %>% 
        extract2(1),
      cont2 = cont_2$cont2 %>% 
        levels %>% 
        extract2(1),
      x = .35,
      y = -.6,
      lab = c("Respondents\nless accurate")
    )
    ) +
  geom_point(
    aes(party_num, estimate),
    data = c2b,
    shape = 21,
    size = 8.5,
    fill = "grey99",
    color = "grey5"
    ) +
  geom_text(
    aes(party_num, estimate, label = lab),
    data = c2b,
    size = 2.5
    ) +
  geom_text(
    aes(x, y, label = lab),
    size = 2.5,
    hjust = 0,
    angle = 90, 
    lineheight = .75,
    data = data.frame(
      dv = cont_2$dv %>% 
        levels %>% 
        extract2(2),
      cont2 = cont_2$cont2 %>% 
        levels %>% 
        extract2(1),
      x = .35,
      y = .025,
      lab = c("Respondents\nfavor regulation")
    )
    ) +
  geom_text(
    aes(x, y, label = lab),
    size = 2.5,
    hjust = 0,
    angle = 90,
    lineheight = .75,
    data = data.frame(
      dv = cont_1$dv %>% 
        levels %>% 
        extract2(2),
      cont2 = cont_1$cont2 %>% 
        levels %>% 
        extract2(1),
      x = .35,
      y = -.2,
      lab = c("Respondents\noppose regs.")
      )
    ) +
  facet_grid(
    dv ~ cont2,
    scales = "free",
    space = "free_x"
  ) +
  scale_x_continuous(
    expand = c(.08, .08),
    breaks = c(1, 3, 5),
    labels = c("Strong\nDem",
               "Independent",
               "Strong\nRep")
  ) +
  labs(
    x = "Partisanship",
    y = ""
  ) 


l_c <- map2(
  list(
    lm_1,
    lm_2
  ) %>% 
    flatten,
  list(
    d1, d1, d2, d2
  ),
  function(i, j)
    # i <- lm_1[[1]]
    # j <- d1
    i %>%
    emmeans(
      pairwise ~ cond|party_num, 
      data = j,
      at = list(party_num = c(1, 3, 5))
    ) %>%
    use_series(contrasts) %>%
    as.data.frame %>% 
    mutate(
      lab = str_c(
        estimate %>% 
          round(2) %>% 
          as.character %>% 
          str_replace("0.", "."),
        p.value %>% 
          stars.pval(),
        " (",
        SE %>% 
          round(2) %>% 
          as.character %>% 
          str_replace("0.", "."),
        ")"
      ) %>% 
        str_replace_all(
          fixed(". "),
          fixed(" ")
        ),
      party_num = party_num %>% 
        mapvalues(
          c(1, 3, 5),
          c("Democrat",
            "Independent",
            "Republican")
        )
    )
)


lc2 <- l_c %>% 
  bind_rows %>% 
  tbl_df %>% 
  mutate(
    mod = 1:4 %>% 
      rep(each = 9)
  ) %>% 
  select(
    mod, everything()
  )

lc3 <- lc2 %>% 
  ddply(
    .(contrast, party_num),
    function(i)
      i$lab %>% 
      as.character
  )


# table 3

stargazer(
  lm_1,
  lm_2, 
  add.lines = list(
    "Auxiliary Quantitities",
    "Misstatement - Correction",
    c("Democrats",
      lc3[7, 3:6] %>% 
        as.character
    ),
    c("Independent",
      lc3[8, 3:6] %>% 
        as.character
    ),
    c("Republican",
      lc3[9, 3:6] %>% 
        as.character
    ),
    "Items - Correction",
    c("Democrats",
      lc3[4, 3:6] %>% 
        as.character
    ),
    c("Independent",
      lc3[5, 3:6] %>% 
        as.character
    ),
    c("Republican",
      lc3[6, 3:6] %>% 
        as.character
    ),
    "Items - Misstatement",
    c("Democrats",
      lc3[1, 3:6] %>% 
        as.character
    ),
    c("Independent",
      lc3[2, 3:6] %>% 
        as.character
    ),
    c("Republican",
      lc3[3, 3:6] %>% 
        as.character
    )
  ),
  covariate.labels = c(
    "Condition: Misstatement",
    "Condiiton: Correction",
    "Partisanship",
    "Cond: Misstatement x Partisanship",
    "Cond: Correction x Partisanship",
    "Constant"), 
  type = "text",
  # type = "latex",
  no.space = T,
  model.numbers = F,
  initial.zero = F,
  keep.stat = c("n", "rsq"),
  digits = 2,
  column.labels = c(
    "May Experiment",
    "June Experiment"
  ),
  column.separate = c(2, 2), 
  dep.var.labels.include = T,
  dep.var.labels   = 
    c("Polar caps size",
      "Enviro regs.",
      "Paris accords flex",
      "Enviro regs")
)

# Table 1

d3 %>% 
  mutate(
    sample = sample %>% 
      factor(
        c("rdd",
          "turk",
          "combined") %>% 
          rev
      )
  ) %>% 
  select(
    resp, study, cond, sample,
    gender, age, ethnicity, educ, income
  ) %>% 
  gather(
    cat, ans, -c(resp, study, cond, sample)
  ) %>% 
  mutate(
    over = "demos"
  ) %>% 
  bind_rows(
    d3 %>% 
      select(
        resp, study, cond, sample,
        trump_approve, partyid
      ) %>% 
      gather(
        cat, ans, -c(resp, study, cond, sample)
      ) %>% 
      mutate(
        over = "atts"
      )
  ) %>% 
  group_by(
    study, cond, sample,  over, cat, ans
  ) %>% 
  tally %>% 
  na.omit %>% 
  mutate(
    prop = n %>% 
      divide_by(
        n %>% 
          sum(na.rm = T)
      ) %>% 
      multiply_by(100)
  ) %>% 
  ungroup %>%  
  mutate(
    cat = cat %>% 
      factor(
        c("partyid", "trump_approve", "age", "educ", "ethnicity", "gender", 
          "income")
      )
  ) %>% 
  unite(
    cv, over:ans, sep = "_"
  ) %>% 
  na.omit %>% 
  mutate(
    cv = cv %>% 
      factor(
        c("atts_partyid_democrat",
          "atts_partyid_independent",
          "atts_partyid_republican", 
          "atts_trump_approve_approve",
          "atts_trump_approve_disapprove", 
          "demos_age_18-39",
          "demos_age_40-59",
          "demos_age_60+",
          "demos_educ_hsd or less",
          "demos_educ_some college",
          "demos_educ_ba +",
          "demos_ethnicity_white", 
          "demos_ethnicity_black",
          "demos_ethnicity_hispanic",
          "demos_ethnicity_asian", 
          "demos_ethnicity_other",
          "demos_gender_female",
          "demos_gender_male",
          "demos_income_<$50k",
          "demos_income_$50-$100k",
          "demos_income_>$100k")
      )
  ) %>% 
  select(
    study, sample, cond:prop
  ) %>% 
  select(-n) %>% 
  mutate(
    prop = prop %>% 
      round
  ) %>% 
  spread(cv, prop)


# figure 4

d5 <- d4 %>% 
  mutate(
    num_trump_approve = trump_approve %>% 
      mapvalues(
        c("approve",
          "disapprove"),
        c(1, 0)
      ) %>% 
      as.character %>% 
      as.numeric,
    num_gender = gender %>% 
      equals("male") %>% 
      as.numeric,
    num_age = age %>% 
      mapvalues(
        d4$age %>%
          unique,
        c(60, 
          mean(c(40, 59)),
          mean(c(18, 39)),
          NA
        )
      ) %>% 
      as.numeric,
    num_partyid = partyid %>% 
      mapvalues(
        c("democrat",
          "independent",
          "republican"),
        1:3
      ) %>% 
      as.character %>% 
      as.numeric,
    num_educ = educ %>% 
      mapvalues(
        c("hsd or less",
          "some college",
          "ba +"),
        1:3
      ) %>% 
      as.numeric,
    num_income = income %>% 
      mapvalues(
        d4$income %>% 
          levels,
        c(50, 75, 100)
      ) %>% 
      as.character %>% 
      as.numeric
  )

t1 <- d5 %>% 
  select(
    resp, cond, study, sample, 
    starts_with("num_")
  ) %>% 
  gather(
    dv, val, starts_with("num_"), na.rm = T 
  )

t2 <- t1 %>% 
  group_by(
    study, sample, dv
  ) %>% 
  by_slice(
    ~lm(
      val ~ cond, data = .
    )
  )

t3 <- pmap_dfr(
  list(
    t2$study,
    t2$sample,
    t2$dv,
    t2$.out
  ),
  function(i, j, k, l)
    
    # i <- t2$study[[1]]
    # j <- t2$sample[[1]]
    # k <- t2$dv[[1]]
    # l <- t2$.out[[1]]
    
    l %>% 
    emmeans("cond") %>%
    contrast("eff") %>%
    tidy %>% 
    mutate(
      study = i,
      sample = j,
      dv = k
    )
) %>% 
  mutate(
    hi = estimate %>% 
      add(
        std.error %>%
          multiply_by(1.96)
      ),
    lo = estimate %>% 
      subtract(
        std.error %>%
          multiply_by(1.96)
      )
  )

t5 <- pmap_dfr(
  list(
    t2$study,
    t2$sample,
    t2$dv,
    t2$.out
  ),
  function(i, j, k, l)
    
    # i <- t2$study[[1]]
    # j <- t2$sample[[1]]
    # k <- t2$dv[[1]]
    # l <- t2$.out[[1]]
    
    l %>% 
    glance %>% 
    select(p.value) %>% 
    # emmeans("cond") %>%
    # contrast("eff") %>%
    mutate(
      study = i,
      sample = j,
      dv = k
    )
) %>% 
  mutate(
    p_lab = p.value %>%
      round(3) %>% 
      as.character %>% 
      str_sub(2)
  )


t3$contrast %<>% 
  mapvalues(
    t3$contrast %>% 
      levels,
    c("Control condition",
      "Correction condition",
      "Misstatement condition")
  ) %>% 
  factor(
    c("Control condition",
      "Misstatement condition",
      "Correction condition")
  )

t3$dv %<>% 
  mapvalues(
    t3$dv %>% 
      unique,
    c("Age",
      "Education",
      "Gender",
      
      "Income",
      "Partisanship",
      "Trump approval")
  )

t5$dv %<>% 
  mapvalues(
    t5$dv %>% 
      unique,
    c("Age",
      "Education",
      "Gender",
      
      "Income",
      "Partisanship",
      "Trump approval")
  )




t5$sample %<>% 
  mapvalues(
    c("rdd",
      "turk"),
    c("RDD",
      "MTurk")
  )


t3$sample %<>% 
  mapvalues(
    c("rdd",
      "turk"),
    c("RDD",
      "MTurk")
  )

t3$dv %<>% 
  fct_reorder(
    t3$p.value
  )

t5$dv %<>% 
  factor(
    t3$dv %>% 
      levels
  )


t3 %>% 
  ggplot() +
  geom_vline(
    xintercept = 0,
    linetype = "dashed"
  ) +
  geom_label(
    aes(
      y = study,
      label = str_c("p = ", p_lab)
    ),
    size = 3,
    x = 0,
    label.size = 0,
    nudge_y = .25,
    fontface = "italic",
    data = t5,
    fill = "grey95"
    ) +
  geom_pointrangeh(
    aes(xmin = lo, xmax = hi, x = estimate, y = study, group = contrast, fill = contrast),
    shape = 21,
    fatten = 3.5,
    size = .25,
    position = position_dodgev(height = .3)
  ) +
  facet_grid(
    sample ~ dv,
    scales = "free_x"
  ) +
  scale_fill_grey(
    end = .96,
    start = .5
  ) +
  scale_y_discrete(
    breaks = c("one", 
               "two"),
    labels = c("1", "2")
  ) +
  labs(
    x = "Conditional difference to overall mean",
    y = "Study",
    fill = ""
  ) 

